PEC1 Análisis de datos ómicos
Máster interuniversitario de Bioestadística y Bioinformática
Información sobre el dataset
Vamos a utilizar el dataset GSE7614. Este dataset corresponde al estudio Phenotypic and transcriptional response to selection for alcohol sensitivity in Drosophila melanogaster.
| Información del dataset | |
|---|---|
| Organismo | Drosophila Melanogaster |
| Número de muestras | 24 (15 individuos por muestra) |
| Grupos | Macho Control (x2 Rep. Biológicas x2 Líneas) |
| Hembra Control (x2 Rep. Biológicas x2 Líneas) | |
| Macho Sensible (x2 Rep. Biológicas x2 Líneas) | |
| Hembra Sensible (x2 Rep. Biológicas x2 Líneas) | |
| Macho Resistente (x2 Rep. Biológicas x2 Líneas) | |
| Hembra Resistente (x2 Rep. Biológicas x2 Líneas) | |
| Número de características | 18952 |
| Plataforma | [Drosophila_2] Affymetrix Drosophila Genome 2.0 Array |
| GPL1322 | |
| Tipo de muestra | in situ oligonucleotide |
Cargado de datos crudos
En este caso hemos descargado todos los archivos .CEL del dataset GSE7614. También hemos consultado a que grupo pertenece cada una de las 24 muestras que contiene el dataset. Hemos preparado también un archivo que contiene toda la información referente a los grupos y las muestras y lo hemos cargado en memoria.
short_name sigue la convención c(ontrol)/r(esistant)/s(ensible) m(ale)/f(emale). Los números separan las diferentes samples y las grupos de los que salieron (ver información sobre el dataset).
Luego hemos generado un ExpressionFeatureSet juntando los archivos .CEL y annotatedDataFrame que contiene el phenoData
Tenemos en total 535824 filas/sondas y 24 muestras.
Control de calidad de los datos crudos
Para realizar un control de la calidad de los datos crudos vamos utilizar la libreria Array Quality Metrics y affyPLM.
Lo haremos de dos maneras diferentes, primero sin transformar los datos con logaritmo de la intensidad y luego haciéndolo.
Sin transformar
Salida Array Quality Metrics sin transformar
En el gráfico A podemos observar un resumen de todas las muestras y el resultado obtenido del control de outliers. Vemos que según el algoritmo las muestras 1 y 17 podrían presentar algún problema por ello las tendremos en especial consideración cuando observemos los gráficos. En este caso no disponemos de la posibilidad de hacerlo, pero podríamos comprobar las imágenes de los microarrays para ver si ha ocurrido alǵun problema durante su procesado.
En el gráfico B observamos una agrupación por distancias de los microarrays, vemos que hay una separación clara según el sexo de la muestra, a excepción de las muestras 1 y 17 que muestran una distancia un poco atípica se encuentran más distantes del resto de muestras de machos.
En el gráfico C observamos los boxplots de las intensidades de cada uno de los arrays. También aparece marcado el array 17 quizá más disperso que el resto, también aparece marcado el 10, sin embargo no parece que sean muestras particularmente aberrantes.
En el gráfico D podemos observar un análisis por componentes principales, veremos debajo uno como el que se ha realizado durante las prácticas de la asignatura, en este caso vemos una separación clara entre los grupos macho y hembra de las muestras lo cual es esperable y deseable.
En el gráfico E, vemos que las muestra siguen un mismo patrón, aunque la escala hace un poco complicada esta observación, probaremos luego con una función de affysimple.
En el gráfico F, el valor esperable es que la mediana sea constante, sin embargo aparece una pequeña joroba en el lado derecho, pero seguramente esté relacionado con una saturación de la intensidades como indica el propio paquete.
En el gráfico G, podemos observar la comparación de cada uno de los arrays respecto a un pseudo-array que tiene el valor mediano a través del resto de arrays. En este caso lo esperable es que los valores se concentren a lo largo del eje M=0. Sin embargo al no estar transformados los valores aparece agrupados en el margen izquierdo haciendo difícil la visualización.
Vamos ahora a ver un PCA más completo
Vemos una agrupación clara por la condición sexo, no tanto por las condiciones experimentales.
Log-Transform
Array Quality Metrics Output
Al estar transformados con el logaritmo la forma de las gráficas cambia ya que hemos estirado la distribución. El significado de los gráficos es el mismo que en el apartado anterior.
El array 17 sigue apareciendo como un outlier, sin embargo en ninguno de los plots aparece particularmente dramático por tanto lo conservaremos.
affyPLM
Vamos ahora a utilizar la librería affyPLM vamos a ajustar un modelo a nivel de sonda y vamos a comprobar las gráficas RLE y NUSE.
En el caso del RLE vamos a encontrar una medida estandarizada de la expresión, que ha de ser similar a través de arrays, y en el caso de NUSE observaremos los residuos normalizados del modelo cuya distribución ha de ser similar a través de arrays.
No parece que ninguno de los arrays se desvíe particularmente en ninguna de las dos gráficas (quizá la c_1_m_1), si que vemos que hay una diferencia clara entre sexos. Podríamos preguntarnos en este caso si se ha cometido algún error y los machos perteneces a un batch y las hembras a otros, es decir hemos generado algún tipo de confound.
Conclusiones de calidad
- Conservamos todos los arrays
- Podría ser útil revisar el 17
- Podría ser útil consultar si ha habido alguna diferencia sistemática entre los arrays de machos y hembras, batch, fecha, operario,…
Normalización
Procedemos ahora a normalizar los datos con el método RMA.
Durante este paso el algoritmo:
- Eliminará el ruido de fondo
- Normalizará la señal entre arrays, nunca entre condiciones experimentales.
- Colapsará todos los valores de cada grupo de sondas en un único valor por gen.
Una vez normalizado vamos a volver a realizar un análisis de calidad de los datos
Array quality metrics output después de la normalización
Array Quality Metrics Output
Los gráficos son los mismos que en apartados anteriores, ahora vemos una separación mucho más clara en el array de distancias. El array número 17 ahora es más similar al resto.
Solo existen dos detalles destacables, uno es la diferencia entre patrones de intensidad en la figura E entre las muestras machos y hembras, y la joroba que aparece en la figura F, cuando la mediana debería de ser plana, quizá sea debida a la diferencia entre sexos.
Conclusiones normalización
Los datos parecen correctos, por tanto no es necesaria ninguna acción.
Filtrado
Vamos a filtrar las probes que no tengan entrada en y las que corresponden a las sondas de control de Affymetrix. Vamos a hacerlo de manera manual utilizando la base de datos correspondiente al chip: Affymetrix Drosophila Genome 2.0 Array, en el paquete ‘drosophila2.db’. Lo que haremos será eliminar todas las features cuyo entrada EntrezID sea NA.
También vamos a eliminar algunas sondas que tienen varios genes asignados, o los genes que tienen varias sondas asignadas. Puede que exista alguna manera más correcta de realizar este paso, pero por sencillez las vamos a eliminar.
De esta manera pasamos de tener 18952 features a 11294. Lo cual nos ayudará, por ejemplo, reduciendo el impacto de la correción por múltiples comparaciones.
Modelo y contrastes
Estamos ante un diseño de 3x2, en este caso vamos a simular que estamos interesados en las comparaciones:
- Main Effect Macho vs. Hembra
- Single Effect Resistente vs Sensible
- Interacción Sexo vs (Resistente vs Sensible)
Vamos a ajustar un modelo lineal con Limma y la siguiente matriz de diseño
Horizontal axis indicates the term of the design matrix. Vertical axis indicate the short name of the sample
Una vez ajustado el modelo calculamos los siguientes contrastes
Anotación
Anotamos los datos con el paquete drosophila2.db.
Resultados
Genes diferencialmente expresados entre condiciones
Las tablas permiten ordenar y filtrar los genes en función de diversos parámetros. Solamente se incluyen los 1000 genes más diferencialmente expresados por cada comparación.
Corrección por multiples comparaciones False Discovery Rate Benajin y Hochsberg CITA!
Macho vs. Hembra
Resistente vs. Sensible
Interacción Sexo y (Resistente vs. Sensible)
Volcano plot
Aquí presentamos un volcano plot con algunos extras sobre el original, con los p-valores ajustados y sin ajustar. Además se incluye una tabla con el porcentaje de genes diferencialmente expresados en cada uno de los contrastes. Se incluye el corte en p<.05 para valores ajustados y sin ajustar.
Volcano plot modificado. Las líneas horizontales indican p-valor=.05. Las líneas verticales indican el un fold-change de 2.
| contrast | Sin corregir p<.05 | Corregido FDR p<.05 |
|---|---|---|
| int_sex_resvssens | 5.2% | 0.13% |
| male_vs_female | 86% | 86% |
| res_vs_sens | 19% | 4.1% |
| a Redondeado a dos decimales |
Diagrama de Venn
Los resultados que aparecen aquí están corregidos al p<.05, con el p-valor corregido a través de múltiples comparaciones con FDR (Bejamini Hochberg) y se han considerado los contrastes de manera separada para esta correción.
| res_vs_sens | male_vs_female | int_sex_resvssens | |
|---|---|---|---|
| Down | 198 | 5825 | 9 |
| NotSig | 10833 | 1619 | 11279 |
| Up | 263 | 3850 | 6 |
Diagrama de Venn que muestra el número de contrastes significativos encada uno de los contrastes y sus intersecciones.
Tabla de anotaciones
Se ha generado (aunque no sea accesible desde el documento, si que está en el repositorio de git) con todos los genes anotados que han resultado signficativos en cada uno de los contrastes.
Archivo: ‘annotations.html’
TODO ANOTAR LOS GENES CON MAYOR FOLD CHANGE
Perfiles de expresión
Análisis de enriquecimiento
- Para estos análisis hemos seleccionado todos aquellos genes cuyo p-valor corregido con FDR es menor a .1.
- Hemos realizado un análisis por contraste
- Hemos utilizado las bases de datos KEGG y GO
- En el caso de KEGG hemos corregido siempre con un p<.1
- En el caso de GO hemos corregido siempre con un p<.1 y un q<.05
- Hemos rebajado ligeramente el umbral de significación para poder mostrar los gráficos y los tipos de análisis disponibles, de otra manera esta sección quedaría vacía.
- En un análisis real mantendríamos el umbral de signifiación fijado a priori, o marcaríamos el análisis como exploratorios.
- Hemos realizado 10000 permutaciones en los GSE
Enlaces a los tipos de gráficos utilizados:
Gene Ontology
Explicación de para que vale la Gene Ontology
Over-representation analysis
Solo cálcula si es por azar que los genes sobre representados aparecen más en un pathway que en otro, con un fisher o un binomial, por azar deberían de aparecer, pero en realidad aparecen.
Macho vs. Hembra
En este caso aparecen muchos genes relacionados con cada ruta y la visualización se hace un poco complicada. Un upsetplot hubiese sido más útil, pero en este caso el paquete daba un error de compatibilidad que no he sido capaz de solucionar.
Resultados de sobrerrepresentación para GO
Resultados de sobrerrepresentación para GO
Resultados de sobrerrepresentación para GO
Resistente vs. Sensible
Resultados de sobrerrepresentación para GO
Resultados de sobrerrepresentación para GO
Resultados de sobrerrepresentación para GO
Interacción Sexo y (Resistente vs. Sensible)
## [1] NAGene Set enrichment analysis
Aquí el approach es un poco distinto porque tenemos en cuenta si se están up o down regulated todos los metabolitos a cañón ahí dentro.
Macho vs. Hembra
Resultados de GSE para GO
Resultados de GSE para GO
Resistente vs. Sensible
## [1] NAInteracción Sexo y (Resistente vs. Sensible)
## [1] NAKEGG
Over-representation analysis
| Contrast | Rutas sobrerrepresentadas |
|---|---|
| Macho vs. Hembra | 2 |
| Resistente vs. Sensible | 0 |
| Interacción Sexo y (Resistente vs. Sensible) | 4 |
Macho vs. Hembra
Resultados de sobrerrepresetanción para KEGG
Resultados de sobrerrepresetanción para KEGG
Resistente vs. Sensible
## [1] NAInteracción Sexo y (Resistente vs. Sensible)
Resultados de sobrerrepresetanción para KEGG
Resultados de sobrerrepresetanción para KEGG
Gene Set enrichment analysis
En los casos en que hemos encontrado rutas significativas, hemos utilizado el paquete pathway para descargar la ruta de la base de datos de KEGG e incluirla en el documento. Así los genes significativos se pueden observar dentro del contexto de la ruta.
| Contrast | Rutas sobrerrepresentadas |
|---|---|
| Macho vs. Hembra | 0 |
| Resistente vs. Sensible | 2 |
| Interacción Sexo y (Resistente vs. Sensible) | 0 |
Macho vs. Hembra
## [1] NAResistente vs. Sensible
## Picking joint bandwidth of 0.926
Resultados de GSE para KEGG
Resultados de GSE para KEGG
Interacción Sexo y (Resistente vs. Sensible)
## [1] NA